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Abstract 

Using Stokes’ formula, a gravimetric geoid has been calculated for Canada. The input data 
are as follows : 15* x 15 ; block averages were used for Canada and the USA and 1 x 1 
block averages and satellite model (GEM— Tl) provided values for the remaining part of the 
Earth. The geoid was calculated at 6 398 points covering the area within the points 
(A is + west) : 

Pt(40,125); pa(75, 184); p 3 (75,10); and p 4 (40,60) 

The computed geoid refers to the GRS1980 and reaches a local minimum of -47.3 metres 
around the western part of Hudson Bay. A contour map of the geoid is shown. 


Introduction 

There is a renewed interest in geoid determination in connection with the use of the Global Positioning 
System (GPS) for geodetic purposes. The demand in this case is for geoidal height differences over rela- 
tively short distances (measured in kilometres), with an error of a few parts per million relative accuracy. 
If the differences in one region are to be connected to those in another region, the need for absolute geoid 
height will demand a similar accuracy to enable one to join the piece-wise information. Gravity coverage 
over the entire surface of the Earth has been assembled for the computation. Computation points over 
continental Canada have been selected at a nearly equidistant interval. The subject of this short note is to 
report on the geoid computation and show the results in the form of a contour map. Further information 
is available in Nagy (1988). 


Theory 


The geoid was calculated from the formula derived by Stokes(l849) : 
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Stokes 9 function is defined as : 
S(ip) = cosec 1 0 


is the geoidal height, 

is the mean radius of the earth, 

is the mean gravity, 

is the gravity anomaly, corresponding to da } 

is the Stokes function, 

is the surface element of the unit sphere. 
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Thus to calculate the geoidal height, N, the gravity anomaly, A g, (representative for the surface 
element), is multiplied by the Stokes function and the area of the surface element, da. This product is 
then summed up over the entire surface of the Earth. This process requires gravity everywhere. In the 
following, the input data used in the computations will be described briefly. 
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Input data 

For the computations, the following data-sets were assembled : 

CANUS15 : 15'x 15' block surface gravity averages for Canada and the continental 
United States, 

WORLD1 : 1° X 1° block surface gravity averages outside the CANUS15 data-set, 

MISSING : 1° x 1° block satellite gravity model values calculated at the centre 
of each block where terrestrial gravity coverage is missing. 

Each of these data-sets will be discussed briefly and some values associated with each particular data- 
set will be given in Table 1. Data-sets not originally on the GEODETIC REFERENCE SYSTEM 1980 
(Geodesist’s Handbook 1980), were transferred to this reference surface durmg the computations. 

CANUS15 

This 15' data-set consists of two files : 

CAN15 : 15' x 15' block surface gravity averages for Canada, 

USA15 : 15' x 15' block surface gravity averages for the continental 
United Statas. 

The CAN15 was calculated from 602 577 point values, available at the begining of 1988 from our 
National Gravity Data Base. The USA15 block averages were obtained from the National Geodetic 
Survey of Rockville, Md. Computed from 1256119 point values, 44029 blocks were made available^ 
Parts of this data-set covering Canada were screened out. Also omitted were some offshore values and 
data over Hawaii. After this data editing, 21 510 blocks were included to represent gravity over the United 
States. The combination of these files resulted in 65 539 15' blocks and provides the best possible surface 
gravity coverage for North America at this resolution. 


WORLD1 

This 1° X 1° mean world-wide surface gravity data-set was obtained from The Ohio State University 
and is described in detail in Despotakis(1986). It contains 48955 values of which 44 203 were used (over 
North America, the more detailed CANUS15 data-set, described earlier, was used). 


MISSING 

After assembling the two previously described files, there were still 15 608 1° x 1° bloc £y™ ] 
surface with no surface gravity values. For these missing blocks the most recent GEM-T1 satellite 
gravity model values were used. The spherical harmonic coefficients were obtained from the Goddard 
Space Flight Center (Marsh et al. 1987) and were evaluated at the center of each block. 

This completed the preparation of the input, consisting of three files, and provided the required gravity 

coverage over the entire surface of the Earth. 

To conclude this section, Table 1 summarizes some statistics about the input files. 


Table 1 : INPUT FILE STATISTICS 



CAN15 

USA15 

CANUS15 

CANlP 

CANlR 

WORLDl 

MISSING 

min 



-223.3 

-112.0 



IHEIsSfil 

max 



239.4 

98.3 




N 



65 539 

3446 

3 446 


n 
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In Table 1, min and max give respectively, the minimum and maximum value of gravity, for the particular 
file and N designates the number of blocks in these files. CAN1P and CAN1R refer to 1 x 1 files 
calculated for Canada : CANlP was obtained simply by averaging the point values within each 1° x 1° 
block, CAN1R is the result of averaging all of the available 15' x 15' mean values within each 1° x 1° 
block with unit weight (up to 16 blocks). This latter value may be more of a regional representation for 
the block than the average of all values in the block. 


Computation 

Early estimates of computing time for the CYBER 730 indicated that 100 CPU hours may be needed to 
compute the geoid at 1 000 points. For a few thousand points, the use of the CYBER was not feasible. 
In the meantime, the Cray 1-S supercomputer became available for use and all the required programs 
were transported to this computer. Further program developments with optimization were then carried 
out in this new environment. For final comparison, the program as used on the Cray was transported 
back to the CYBER and run for performance comparisons. The test computations indicated a 34-fold 
speed increase on the Cray. 

The computations were always carried out in such a manner that, for each computation point the 
contribution of each of the three data-sets was also available separately. This allowed each effect to be 
seen separately when plotted as a partial contribution to the geoidal height. The computation points 
were selected first at 1° intervals along the meridian and at the same distance along the parallels. For a 
denser net of computing points, this grid was shifted in both directions. Additional points were inserted 
where better definition of the geoid was required. Finally about 1000 points were added to trace out 
the geoid locally in greater detail. This resulted in a total of 6 398 computation points, from which the 
present geoid map was prepared. In addition to a colour APPLICON map, a contour map was also 
prepared, and shown in Figure 1. As mentioned earlier the geoid refers to the GRS1980. It reaches a 
local minimum of -47.3 m at <p = 59.76° and A = 92.34° W in Hudson Bay. 

The total computation time on the Cray wag 8.07 CPU hours, which corresponds to over 275 hours 
on the CYBER. 
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